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ABSTRACT 

Initialization  of  a  five-level  global  baroclinic 
primitive  equation  model  was  examined  using  real  data. 
Experiments  were  conducted  using  a  Robert  (1965)  time 
frequency  filter,  Euler  backward  time  integration,  and  an 
iterative  initialization  scheme  to  determine  the  effect  of 
each  on  the  generation  of  inertial-gravity  waves  resulting 
from  an  improper  balance  between  initial  mass  and  wind 
fields.   In  addition,  a  global  sigma-surface  linear 
balance  equation  solution  was  examined  for  the  purpose  of 
determining  initial  winds.   Although  all  forecast  fields 
were  meteorological,  certain  undesirable  features  were 
generated  in  the  polar  and  equatorial  regions. 
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I.   INTRODUCTION 

The  proposed  FNWC  five-level,  global,  weather  prediction 
model  using  a  spherical,  staggered,  sigma  coordinate  system 
is  currently  under  development  at  FNWC  to  augment  or  replace 
FNWC  Northern  Hemispheric  model  on  an  operational  basis. 
The  model  was  converted  to  operate  on  the  IBM  360  at  the 
W.  R.  Church  Computer  Center  for  continued  student  and 
faculty  research  at  the  Naval  Postgraduate  School. 

The  purpose  of  this  particular  study  was  to  examine 
various  schemes  for  initialization  of  the  model  using  real 
data  from  FNWC  objective  analyses.   A  previous  study  by 
Elias  (1973)  and  a  current  study  by  Mihok  (1974)  utilized 
analytically  derived  fields  to  examine  the  model  charac- 
teristics.  The  use  of  real  data  permitted  real,  operational 
type  problems  to  be  present  and  examined  in  the  model 
initialization. 

The  objective  was  to  determine  the  best  method  of 
controlling  inertial-gravity  waves  generated  by  the  initial 
imbalance  between  the  mass  and  wind  fields  with  the  aim  of 
improving  the  model  in  the  early  stages  of  a  forecast. 
These  spurious  disturbances  normally  diminish  due  to  the 
"geostrophic  adjustment"  mechanism  inherent  in  the  model  so 
that  the  24-hour  forecasts  and  longer  are  not  seriously 
affected.   However,  as  satellite  data,  which  are  essentially 
continuous  in  time,  are  increasingly  used  on  a  regular  basis 
between  synoptic  times  to  update  numerical  predictions,  it 
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will  be  highly  desirable  to  suppress  this  inertial-gravity 
"noise"  resulting  from  the  introduction  of  new  data. 

A  procedure  was  developed  to  obtain  initial  balanced 
wind  fields  on  sigma  surfaces  for  comparison  with  the 
balanced  winds  on  pressure  surfaces. 

In  addition,  tests  were  made  with  various  combinations 
of  the  Robert  time  frequency  filter,  periodic  usage  of  an 
Euler  backward  time  step,  the  sigma-surface  balanced  winds, 
and  a  dynamic  balancing  scheme  suggested  by  Temperton  (1973) 
in  order  to  determine  the  best  method  of  controlling  the 
initial  trauma. 
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II.   BAROCLINIC  PRIMITIVE  EQUATION  MODEL 

The  governing  equations  in  the  model  (Elias,  1973)  are 
given  below.   These  equations  are  similar  to  those  used  by 
Smagorinsky,  et  al .  (1965),  Arakawa,et  al.  (1969)  and 
Kesel  and  Winninghoff  (1972).   A  brief  description  of  the 
physics  of  the  model  may  be  found  in  detailed  form  in  Kesel 
and  Winninghoff  (1972).   Mihok  (1974)  lists  a  complete  set 
of  finite  difference  equations  in  spherical  form  used  in 
the  current  model.   Appendix  A  contains  the  finite  differ- 
ence forms  of  the  sigma  coordinate  linear  balance  equation, 
the  vector  form  is  given  below. 

A.   PRIMITIVE  EQUATIONS 

The  vector  form  of  the  equations  of  motion  in  the  sigma 
coordinate  system  are  as  follows: 

Momentum  equation 

ttt(ttV)    +    (v-TrV)V  +    TT^-(aV)    +    f(KXTrV)    +   ttV<{>   +   ottciVtt 

-  F   +   D 
m 

Thermodynamic   equation 


^(ttT)    +    V-(ttTV)    +    TT^r(Ta)    -   ~(a~  +    cWtt   +    tic)    =    itH   +   Dt 

P 
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Mass  continuity  equation 


|f+  V.(wV)  t.Jf-0 


Moisture  continuity  equation 


^r(irq)  +  V-(irqV)  +  ir~(qa)  =  ttQ  +  D 


Hydrostatic  equation 


u +  ™  -  ° 

3a 


Equation  of  state 

aP  =  RT 

The  vertical  coordinate  sigma  (a)  is  related  to  pressure 
(P)  and  surface  pressure  (tt)  by  the  equation 

0  =  P/tt 

The  terms  and  symbols  are  identical  to  those  used  in 
Kesel  and  Winninghoff  (1972).   A  List  of  Symbols  and 
Abbreviations  appears  on  page  8  of  this  paper. 

B.   GRID 

The  grid  used  by  the  model  is  a  staggered,  spherical, 
five-degree,  latitude-longitude  grid.   Mass  variables 
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(<}),  T,  q,  w)  are  carried  at  10-degree  intervals  over  the 
globe  and  at  the  poles.   Motion  variables  (u,v)  are  carried 
at  alternate  10-degree  intervals  over  the  globe.   Figure  1 
depicts  the  grid  configuration. 

C.  VERTICAL  LAYERS 

Figure  2  depicts  the  five  layers  vertically  distributed 
in  the  model  and  the  variables  associated  with  each  layer. 
Specific  humidity  (q)  is  carried  in  only  the  three  lowest 
layers,  at  the  .9,  .7,  and  .5  sigma  levels.   The  remaining 
variables  u,  v,  T,  and  <J>  are  carried  at  the  .9,  .7,  .5 ,    .3, 
and  .1  sigma  levels.   The  vertical  motion  variable,  w,  is 
calculated  diagnostically  from  the  continuity  equation  at 
layer  interfaces,  the  .8,  .6,  .4,  and  .2  sigma  levels. 

D.  TIME  DIFFERENCING 

The  two  basic  difference  schemes  followed  are  the 
centered  difference  scheme  (leapfrog)  and  the  Euler  backward 
difference  scheme.   The  centered  difference  form  is 


F(t  +  At)  =  F(t  -  At)  +  2At9^t) 


Since  this  difference  form  is  of  higher  order  than  that  of 
the  differential  equation  a  computational  mode  exists  and 
solution  separation  at  even  and  odd  time  steps  will  develop, 
perhaps  even  leading  to  a  type  of  Instability  (Haltiner,  1971) 
In  addition,  the  leapfrog  method  cannot  be  used  on  the  first 
time  step. 
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To  eliminate  the  solution  separation  and  also  to  initi- 
ate the  integration  procedure  the  Euler  backward  scheme  is 
used: 


F*  =  F(t)  +  At  -P(t) 


3t 


F(t  +  At)  =  F(t)  +  At  |~ 


The  Euler  backward  scheme  is  applied  at  six-hour  intervals 
in  most  of  the  experiments  conducted,  however  where  indi- 
cated the  scheme  was  applied  more  often. 

A  time  step  of  ten  minutes  was  used  in  all  experiments. 
In  order  to  insure  linear  computational  stability  (Haltiner, 
1971) >  a  much  shorter  time  step  would  have  to  be  taken  in 
a  latitude-longitude  grid  where  for  a  fixed  longitude 
increment  the  east-west  distance  between  grid  points  decreases 
toward  the  poles.   To  avoid  the  shorter  time  step  the  model 
employs  the  Arakawa  (Gates,  1971)  technique  of  filtering 
or  damping  short  waves  that  would  lead  to  instability  near 
the  poles  thus  allowing  the  10-minute  time  step  throughout. 
Specifically  a  weighted  average  is  taken  of  those  longitu- 
dinal derivatives  that  are  involved  in  gravity  wave  propaga- 
tion with  increasing  smoothing  toward  the  poles.   This  tech- 
nique appears  to  be  superior  to  decreasing  the  time  step 
or  skipping  grid  points  toward  the  poles. 

The  heating  package  is  called  every  six  time  steps  while 
lateral  diffusion,  friction,  convectlve  adjustment,  convective 
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condensation,  and  large  scale  condensation  are  computed 
every  time  step. 
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III.   INITIALIZATION 

A.  DATA  INPUT 

The  real  data  used  in  the  experiments  was  provided  on 
magnetic  tape  by  FNWC  from  objective  analyses  on  the 
Northern  Hemisphere  63  x  63  point  grid. 

The  model  interpolates  the  63  x  63  data  to  a  latitude- 
longitude  grid  of  five-degree  intersections  in  the  Northern 
Hemisphere.   The  grid  interpolation  values  are  then  reflected 
in  the  Southern  Hemisphere  to  allow  the  model  to  integrate 
over  the  globe.   FNWC  is  currently  converting  the  model  to 
read  data  directly  on  the  global  latitude-longitude  grid 
with  real  input  to  the  Southern  Hemisphere. 

In  all  tests  data  consisted  of  12  constant-pressure 
level  temperature  analyses,  10  height  analyses,  4  moisture 
analyses,  sea-level  pressure  and  sea-surface  temperature 
analyses.   The  date  of  the  data  for  the  experiments  is 
12  Z  10  May  1972. 

The  terrain  height  provided  to  the  model  is  an  unsmoothed 
"gross"  terrain.   In  addition  the  interpolation  to  the 
latitude-longitude  grid  yields  an  unrealistic  height  for 
the  interpolated  grid  point.   The  conversion  to  global  input 
should  eliminate  this  problem. 

B.  LINEAR  BALANCE  ON  PRESSURE  SURFACES 

The  stream  function  field  is  derived  from  the  solution 
to  the  linear  balance  equation 


18 


V2^  =  -  Iv^'Vf  +  jrV2<f>  (1) 


using  height  values  horizontally  interpolated  from  the  polar 
stereographic  63  x  63  grid  to  the  spherical  grid.   Equation 
(1)  was  solved  using  the  over-relaxation,  iterative  tech- 
nique (Haltiner,  197D  with  a  modified  coriolis  parameter 
in  the  equatorial  region.   The  wind  fields  were  then  derived 
from  the  computed  stream  function  field  as  follows: 


u  =  -  i  $£  (2) 

u     a  36  Kd) 


84  (3) 


a  cos  6  3A 

After  the  u  and  v  components  have  been  computed  on  the 
pressure  surface  spherical  coordinate  system,  the  components 
are  then  interpolated  vertically  to  the  sigma  coordinate 
system  through  a  logarithmic  interpolation  scheme.   The  u 
and  v  wind  fields  are  then  reflected  to  the  Southern 
Hemisphere . 

C.   LINEAR  BALANCE  ON  SIGMA  SURFACES 

Sundquist  (1973)  stated  that  the  interpolation  process 
described  in  the  previous  section  by  the  very  process  of 
interpolation  further  stimulates  the  inertial-gravity 
oscillations  in  the  early  stages  of  the  integration.   By 
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balancing  the  mass  and  stream  function  fields  directly  on 
the  model  sigma  surfaces  inertial-gravity  waves  were 
suppressed  (Sundquist,  1973).   The  non-linear  balance 
equation  (Sundquist,  1973)  is: 


U  ^     l3xL3y  9x  Kv    dyJ  9x  gy  Kv    gy;j 


3y    L8x    3y    ^u    3x;         9y    3x    ^    3x;j 


=   V-[ttVcJ>  +   RTVtt]  (4) 

Due  to  time  limitations  the  complete  balance  equation  was 
not  tested;  however  the  following  linear  form  was  solved 
over  the  entire  globe: 


fV2i|;   +    (Vip-Vf)    =   ttV24>   +   Vtt-V<}>  +   R[TV2tt   +   VT-Vtt]  (5) 


The  finite  difference  form  of  this  equation  is  given  in 
Appendix  A.   The  Fortran  program  used  to  solve  Equation  (5) 
is  on  page  85  . 

The  correct  value  of  the  coriolis  parameter  was  used 
everywhere  except  in  the  equatorial  region  where  two  alterna- 
tives were  tested  to  obtain  convergence  of  the  relaxation 
scheme.   The  initial  guess  for  the  \\>   field  was  made  with 
the  coriolis  parameter  held  constant  from  10  degrees,  in 
one  case  and  20  degrees  In  the  other  case,  to  the  equator 
in  both  hemispheres.   Elsewhere  the  correct  value  of  the 
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coriolis  parameter  was  used  at  both  north  and  south  latitudes, 
initially  and  during  relaxation.   This  scheme  yielded  a 
positive  stream  function  in  Northern  Hemisphere  and  a  nega- 
tive function  in  the  Southern  Hemisphere.   During  relaxa- 
tion the  value  of  ip  at  the  equator  was  obtained  by  averaging 
the  values  at  5N  and  5S. 

D.  TIME  FREQUENCY  FILTER 

The  Robert  (1965)  time-frequency  filter  incorporated 
into  the  model  utilizes  the  following  averaging  operator: 

F(t)*  =  F(t)  +  a[F(t-l)  +  F(t+1)  -  2F(t)]     (6) 

where  a  is  the  filter  parameter,  t  is  time  and  F  is  the 
field  being  filtered.   Different  values  of  a  were  tested 
ranging  from  0  (no  filtering)  to  .5  (maximum  filtering). 
This  filter  tends  to  damp  the  inertial-gravity  waves 
arising  from  the  imbalance  between  the  mass  and. wind  fields, 
especially  2At  periods,  and  also  essentially  eliminates 
solution  decoupling  (Robert,  1965). 

E.  ITERATIVE  AVERAGING  OF  WIND  COMPONENTS 

An  iterative  scheme  of  averaging  the  u  and  v  components 
of  the  wind  field  derived  from  both  a  forward  and  a  backward 
integration  of  the  model  for  specified  periods  of  time  such 
as  1,  3>  and  6  hours  was  developed.   The  components  derived 
from  each  set  of  forward  and  backward  Integrations  were 
separately  averaged  to  constitute  new  initial  winds  to  begin 
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the  Integration  from  time  zero.   On  the  other  hand,  the 
mass  and  temperature  fields  were  restored  to  their  original 
values  at  time  zero.   This  procedure  was  repeated  two  to 
six  times. 

The  backward  integration  used  a  ten  minute  time  step. 
Due  to  the  non-reversibility  of  the  heating  routine,  dia- 
batic  heating  was  not  incorporated  in  the  backward  integra- 
tion. 

Winninghoff  (1971)  and  Temperton  (1973)  have  shown 
this  type  of  procedure  is  quite  effective  in  removing  the 
imbalance  between  the  initial  mass  and  wind  fields.   The 
process,  as  they  applied  it,  is  too  time  consumming,  however, 
when  considered  from  an  operational  point  of  view. 
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IV.   RESULTS 

Eight  experiments  will  be  presented  in  this  section. 
The  first  three  experiments  were  performed  to  determine 
what  filter  parameter  value  would  be  used  in  subsequent 
experiments  and  how  often  the  Euler  backward  time  step 
would  be  used.   Experiment  IV  consisted  of  two  12-hour 
forecasts  using  an  iteration  initialization  scheme. 
Experiments  V  and  VI  involved  a  comparison  of  winds  derived 
from  a  sigma-surface  balance  equation  and- winds  derived 
from  a  pressure-surface  balance  equation  both  without 
terrain.   Experiment  VII  introduced  terrain  into  the  sigma- 
surface  balance  scheme.   Experiment  VIII  consisted  of  an 
84-hour  forecast  using  a  filter  parameter  of  0.4. 

In  certain  indicated  experiments  plots  of  terrain 
pressure  versus  forecast  time  for  12  hours  were  made  at 
three  latitudes,  20  N,  40  N,  and  80  N  on  110  E  longitude. 
The  small  amplitude  waves  depicted  are  a  reflection  of  the 
inertial-gravity  waves  generated  by  the  lack  of  a  proper 
balance  between  the  mass  and  wind  fields.   As  a  consequence 
of  the  "geostrophic  adjustment"  mechanism  inherent  in  the 
model  the  waves  gradually  tend  to  diminish.   The  degree  of 
time  smoothing  has  a  noticeable  effect  on  these  initial 
oscillations.   In  addition,  the  iteration  initialization 
scheme  also  has  a  damping  effect  on  the  initial  oscillations, 
particularly  on  the  longer  periods.   Although  the  initial 
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winds  obtained  from  the  sigma-surface  linear  balance  equa- 
tion were  expected  to  help  damp  initial  inertial-gravity 
wave  oscillations,  the  results  do  not  demonstrate  the 
damping. 

Chart  A  depicts  the  surface  pressure  objective  analysis 
from  PNWC  magnetic  tape  used  as  initial  input  for  the  model. 
Chart  B  is  also  of  the  initial  surface  pressure  but  is  a 
result  of  the  model  interpolating  from  the  Northern  Hemi- 
sphere 63  x  63  grid  to  the  latitude-longitude  grid  and 
back  to  the  63  x  63.   Some  original  detail  in  the  surface 
pressure  field  is  lost  in  the  process. 

1.   Experiment  I 

Four  12-hour  forecasts  were  made,  each  forecast 
using  a  different  filter  parameter.   The  Euler  backward  time 
step  was  used  once  every  six  hours  in  each  case.   The  fore- 
cast surface  pressure  fields  resulting  from  use  of  filter 
parameter  values  of  0.0,  0.3,  0.4,  and  0.5  can  be  found  on 
Charts  C-F,  respectively.   Figures  3-5  show  four  curves  per 
figure,  each  curve  a  plot  of  terrain  pressure  from  each  of 
the  forecasts  versus  forecast  time  for  the  three  latitude 
grid  points  previously  mentioned.   Curve  1  in  each  figure 
demonstrates,  the  unfiltered  case  and,  as  a  consequence,  the 
most  oscillations.   Filtering  very  effectively  damps  short 
period  oscillations,  as  is  evident  from  Curves  2,  3,  and  4 
in  each  figure.   In  Figure  5,  the  terrain  pressure  at  80N 
tends  to  increase  as  the  filtering  is  increased.   A  value 
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of  a  -  0.4  was  selected  as  a  reasonable  compromise  with 
respect  to  limiting  the  rise  of  pressure  at  the  poles  and 
still  accomplishing  a  significant  damping  of  the  inertial- 
gravity  noise. 

2.  Experiment  II 

Three  12-hour  forecasts  were  produced  in  this 
experiment  using  .the  Euler  backward  time  step  at  1,  2,  and 
3  hour  intervals.   The  forecast  surface  pressure  fields  are 
depicted  on  Charts  G-I.   In  addition  Chart  C  depicts  the 
surface-pressure  forecast  using  the  Euler  backward  time 
step  every  six  hours.   In  this  experiment  the  filter 
parameter  a  was  set  to  zero  so  that  no  filtering  was 
performed.   Figures  6-8  show  the  terrain  pressure  versus 
forecast  time  for  each  forecast  at  the  three  latitudes 
selected.   Curve  1  in  each  figure  represents  the  case  where 
a  =  0.0  and  the  Euler  backward  is  used  once  every  six 
hours.   Curves  2,  3,    and  4  represent  the  Euler  step  being 
used  every  one,  two,  and  three  hours,  respectively. 
Although  very  short  period  oscillations  are  damped,  the 
longer  period  oscillations  that  were  significantly  damped 
by  the  time  frequency  filter  are  not  appreciably  affected. 

3.  Experiment  III 

Combinations  of  a  filter  value  a  =  0.4  and  the 
Euler  backward  time  step  every  one  hour  and  every  six  hours 
were  used  to  produce  12-hour  forecasts.   Figures  9-11 
depict  three  curves  each  of  terrain  pressure  versus  forecast 
time;  Curve  1  representing  a  =  0.4  and  1  hour  use  of  the 
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Euler  step  while  Curve  2  represents  0=0.4  and  6  hour  use 
of  the  step.   Curve  3  was  added  for  comparison  and  repre- 
sents a  =  0.0  or  no  filtering  and  the  Euler  backward  time 
step  used  once  every  six  hours.   As  can  be  readily  seen  at 
all  three  latitudes,  increasing  the  use  of  the  Euler 
backward  time  step  does  not  appreciably  alter  the  degree 
of  damping  of  oscillations  when  used  in  conjunction  with 
the  filter  parameter  a  =  0.4. 
4 .   Experiment  IV 

Two  12-hour  forecasts  were  made  utilizing  an  itera- 
tion initialization  equivalent  in  model  time  to  a  12-hour 
forecast.   Chart  J  depicts  the  surface  pressure  forecast 
obtained  using  a  one-hour  forward/one-hour  backward 
integration  followed  by  averaging  and  repeated  six  times. 
Chart  K  demonstrates  the  surface  pressure  forecast  obtained 
using  the  iteration  scheme  repeated  twice  for  a  three  hour 
forward  and  three  hour  backward  iteration.   Both  forecast 
procedures  used  a  filter  parameter  of  a  =  0.4  and  the 
Euler  backward  time  step  once  every  six  hours.   Figures 
12-14  show  plots  of  terrain  pressure  versus  forecast  time, 
Curve  1  representing  the  one-hour  scheme  and  Curve  2  the 
three-hour  scheme.   Curve  3  represents  the  plot  obtained 
by  using  a  normal  non-iterative  initialization  with  a  =  0.4 
and  the  Euler  backward  step  once  every  six  hours.   Both 
iteration  schemes  damp  the  longer  period  oscillations 
remaining  after  time  filtering. 
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5.   Experiment  V 

Initial  winds  derived  from  the  linear  balance 
equation  solved  on  sigma  surface  without  terrain  were 
compared  with  initial  winds  computed  on  pressure  surfaces. 
The  initial  winds  obtained  directly  on  sigma  surfaces 
were  used  to  integrate  the  model  to  seventy-nine  hours. 
The  12-hour  and  72-hour  forecast  surface  pressure  fields 
are  shown  on  Charts  X  and  Z  to  demonstrate  the  meteorological 
product.   Chart  Y  depicts  the  12-hour  surface  pressure  fore- 
cast field  obtained  using  initial  winds  derived  from-  the 
linear  balance  equation  solved  on  pressure  surfaces  and 
interpolated  to  sigma  surfaces.   Again  terrain  was  not  used 
in  order  to  hopefully  show  similarity  between  the  stream 
function  fields  produced. 

Charts  L  and  M  show  the  initial  stream  function 
fields  on  the  0.9-sigma  level  and  925-mb  level,  respectively. 
Charts  N  and  0  show  the  final  iji  fields  for  the  two  levels. 
Although  the  initial  guess  i|>  fields  are  considerably 
different,  the  final  fields  show  quite  similar  trough  and 
ridge  patterns.   Charts  P  and  Q  show  the  isotach  patterns 
for  the  0.9-sigma  level  and  925-mb  level,  respectively. 
The  magnitude  of  the  winds  on  the  sigma  surface  are  5  to  10 
meters  per  second  slower  than  those  on  the  pressure  surface. 

Charts  R-W  show  the  same  fields  for  the  0.5-sigma 
level  and  500-mb  level  that  Charts  L-P  did  for  the  0.9-sigma 
level  and  925-mb  level.   Again  the  final  ifi  fields  have 
similar  trough  and  ridge  patterns  although  the  first  guess 
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fields  are  different.   Also  the  wind  intensity  on  the 
0.5-sigma  surface  is  less  than  that  on  the  500-mb  surface. 

Figures  15-17  show  two  curves  of  terrain  pressure 
versus  forecast  time  for  the  two  balance  systems.   Curve  1 
represents  the  pressure  case  while  Curve  2  represents  the 
sigma  case.   In  both  cases  a  =  0.4  was  used  with  an  Euler 
backward  time  step  once  every  six  hours. 

In  this  particular  test  the  initial  guess  for  the 
if)  field  in  the  sigma-surface  balance  program  was  made  with 
a  constant  coriolis  parameter  from  10  degrees  to  the 
equator  in  both  hemispheres.   The  correct  coriolis  value 
was  used  everywhere  else,  initially  and  during  the 
relaxation. 

At  the  equator  where  f  =  0,  the  balance  equation 
is  no  longer  of  the  Poisson  or  Helmholtz  type  to  which  the 
relaxation  scheme  is  applicable.   To  avoid  the  difficulty, 
the  value  of  tJ;  at  the  equator  was  obtained  by  averaging 
the  value  at  5N  and  5S  during  each  pass  in  the  relaxation 
scheme . 

6 .   Experiment  VI 

Another  initial  guess  for  the  stream  function  in 
the  sigma-surface  balance  scheme  was  made  in  this  experiment 
The  value  of  the  coriolis  parameters  was  held  constant  this 
time  from  20  degrees  to  the  equator  in  both  hemispheres. 

Again,  the  correct  coriolis  parameter  was 

used  everywhere  else  initially  and  during  relaxation. 
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The  equator  value  of  i|>  was  again  obtained  by  averaging  the 
values  at  5N  and  5S  during  the  relaxation. 

Chart  AA  depicts  the  .9-sigma  level  final  stream 
function  field.   Chart  BB  shows  the  12-hour  forecast 
surface  pressure  field.   Although  Chart  AA  and  N  should 
be  identical  if  convergence  to  the  same  solution  was  being 
accomplished,  it  is  evident  from  Charts  BB  and  X  that  the 
solution  to  the  linear  balance  equation  produced  is  suffi- 
cient to  yield  similar  surface  pressure  forecasts. 

7.  Experiment  VII 

This  experiment  is  identical  in  conditions  to  the 
sigma-surface  balance  part  of  Experiment  V  except  that 
terrain  was  included.   As  previously  stated,  the  terrain 
currently  used  in  the  model  is  an  unsmoothed  terrain. 
Although  the  introduction  of  terrain  to  a  scheme  that  is 
intended  to  take  advantage  of  the  terrain  in  balancing  the 
mass  and  wind  fields  is  probably  not  the  sole  cause,  the 
model  immediately  went  unstable  in  3  hours  and  a  forecast 
was  not  obtained. 

8 .  Experiment  VIII 

An  84-hour  forecast  was  made  with  the  model  using 
a  filter  parameter  a  =  0.4  and  the  Euler  backward  time  step 
once  every  six  hours.   Charts  CC  and  EE  show  the  72-hour 
and  84-hour  forecast  surface  pressure  fields  and  charts  DD 
and  FF  show  the  FNWC  verifying  analyses.   Of  note  are  the 
extreme  pressure  gradients  generated  in  the  polar  region. 
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VI.   CONCLUSIONS 

Various  schemes  for  initializing  a  five-level,  global, 
weather  prediction  model  using  a  spherical,  staggered, 
sigma  coordinate  system  were  examined  in  order  to  determine 
an  optimum  method  for  controlling  inertial-gravity  waves 
generated  by  the  initial  imbalance  between  the  mass  and 
wind  fields. 

The  Robert  (1965)  time  frequency  filter  was  quite 
successful  in  damping  the  inertial-gravity  waves,  expecially 
the  computational  mode  which  has  a  period  of  approximately 
2At.   Values  of  a  from  0.0  to  0.5  were  tried  in  early 
experiments.   The  value  of  a  =  0.4  appeared  the  best 
compromise  between  rising  polar  terrain  pressure  and 
damping  and  was  used  in  subsequent  experiments.   Since  the 
time  frequency  filter  essentially  eliminates  solution 
decoupling,  the  need  for  restarting  with  the  Euler  backward 
time  step  is  eliminated.   It  would  still  serve  a  useful 
purpose  for  the  first  one  or  two  time  steps.   More  frequent 
use  does  not  seem  warranted  since  damping  due  to  use  of  the 
Euler  scheme  is  minimal  compared  to  the  time  frequency 
filter. 

The  iterative  scheme  (Temperton,  1973)  of  averaging  the 
u  and  v  components  of  the  wind  fields  derived  from  a 
forward  and  backward  integration  of  the  model  further 
damped  the  longer  period  oscillations  already  partially 
damped  by  the  Robert  filter. 
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Additional  work  and  testing  is  required  on  the  solution 
of  the  balance  equation  on  sigma  surfaces.   Although 
convergence  to  the  program  criteria  was  indicated,  simply 
changing  the  initial  guess  for  the  stream  function  field 
altered  considerably  the  final  field,  implying  that  true 
convergence  toward  the  correct  solution  was  not  accomplished 
by  the  program.   Nevertheless  the  similarity  between  the 
final  stream  function  fields  produced  on  pressure  surfaces 
and  those  produced  without  terrain  on  sigma  surfaces 
warrants  further  program  testing.   Moreover,  actual  con- 
vergence of  the  iterative  procedure  to  a  correct  solution 
of  the  balance  equation  may  not  be  necessary  or  even  desir- 
able if  several  iterations  will  produce  better  initial  winds 
than  the  pressure-surface  balanced  winds. 

The  inclusion  of  terrain  in  the  sigma-surface  balance 
scheme  caused  instability  in  the  model  very  rapidly.   As 
previously  mentioned,  the  terrain  used  in  the  model  is  not 
smooth.   Smoothing  may  reduce  the  trauma  encountered  with 
the  sigma-surface  balance  program  integration,  but  it  is 
not  conclusive  at  this  time  if  terrain  is  the  sole  cause 
of  instability. 

Perhaps  the  solution  of  the  full  balance  equation  will 
reduce  the  inertial-gravity  waves  at  the  start  of  the 
integration  as  indicated  by  Sundquist  (1973) »    although  the 
linear  balance  version  was  not  effective  in  these 
experiments . 
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The  problems  encountered  with  the  linear  balance  equation 
in  the  equatorial  regions  where  the  coriolis  parameter  Is 
negligible  and  pressure  gradients  are  small  could  be 
resolved  through  another  approach  to  the  problem.*  Using 
the  FNWC  objectively  analyzed  winds  for  the  equatorial 
region  for  instance  from  40S  to  40N  the  vorticity  at  each 
grid  point  could  be  calculated  using: 


3v  .     8u  .     v  . 
_   obs     obs  ,   obs  ,    Q  ,„x 

?obs  -  — 93E 37—  +  "a"  tan  6  (7) 


The  stream  function  ty   could  then  be  determined  by  solving 
the  Poisson  equation: 


V'*  '  <=obs  <8) 


The  calculated  values  of  ty ,    say  at  20  or  30  N  and  S  could 
then  be  used  as  boundary  conditions  to  solve  the  non-linear 
balance  equation  for  the  stream  function  in  the  regions 
poleward  of  the  selected  latitudes  from  the  analyzed 
geopotential  fields  in  the  standard  manner.   Continuity  of 
the  stream  function  between  the  equatorial  belt  and  the 
regions  poleward  would  be  automatically  satisfied.   Finally, 
the  mass  field  (4>)  in  the  tropical  belt  should  be  obtained 
by  solving  the  balance  equation  using  the  previously 
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calculated  stream  function  together  with  the  objectively 
analyzed  <J>  values  at  the  boundaries.   The  u  and  v  irrota- 
tional  wind  components  are  then  derived  from  the  final 
global  stream  function  field. 

The  long  term  integration  indicates  a  need  for  a 
review  of  the  polar  region  finite  differencing  and  the 
method  of  handling  the  pole  points.   Severe  "noise" 
develops  at  the  pole  causing  model  instability  at  about 
88  hours  prognostic  time. 

Although  it  is  evident  that  additional  experimentation 
is  desired  for  the  purpose  of  removing  initialization 
trauma,  the  reader  should  bear  in  mind  that  in  these 
experiments  the  method  of  data  input  by  interpolating  to 
latitude-longitude  grid  points  from  a  polar  stereographic 
grid  probably  introduces  further  imbalance  between  mass  and 
wind  fields. 

On  the  basis  of  these  experiments  a  definite  conclusion 
reached  is  that  for  the  purpose  of  reducing  inertial-gravity 
wave  noise,  a  combination  of  the  Robert  time  frequency 
filter  (with  a  =  0.4)  and  a  forward  and  backward  dynamic 
balancing  scheme  is  quite  effective  and  not  overly  time- 
consuming.   Conversion  to  a  global  latitude-longitude  grid 
data  input  system  coupled  with  a  full  balance  equation 
solved  on  sigma-surfaces  may  perhaps  reduce  the  inertial 
trauma  sufficiently  to  obviate  the  need  for  dynamic  balancing 
by  the  Temperton  (1973)  method. 
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FIGURE  1.   LOCATION  OF  VARIABLES. 

i  is  the  column  counter  and  j  is  the  row  counter. 
Point  (1,1)  is  at  85S,  10E. 
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FIGURE  2.   VERTICAL  LAYERING 

In  the  above  figure,  sigma  (a)  is  the  dimensionless  vertical 
coordinate,  u  and  v  are  the  zonal  and  meridional  wind 
components,  respectively,  q  is  the  specific  humidity,  <J>  is 
the  geopotential,  tt  is  the  terrain  pressure  and  w  is  the 
vertical  velocity  (-a). 
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Minute  at  80M  HOE  (Experiment  I) 
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FIGURE  6.   Plot  of  Terrain  Pressure  Versus  Forecast 
Minute  at  20N  HOE  (Experiment  II) 
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FIGURE  7.   Plot  of  Terrain  Pressure  Versus  Forecast 
Minute  at  HOU   HOE  (Experiment  II) 
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FIGURE  8.   Plot  of  Terrain  Pressure  Versus  Forecast 
Minute  at  80M  HOE  (Experiment  II) 
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FIGURE  9.   Plot  of  Terrain  Pressure  Versus  Forecast 
Minute  at  20N  HOE  (Experiment  III). 
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FIGURE  10.   Plot  of  Terrain  Pressure  Versus  Forecast 
Minute  at  40N  HOE  (Experiment  III) 
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FIGURE  11.   Plot  of  Terrain  Pressure  Versus  Forecast 
Minute  at  80N  110E  (Experiment  III) 
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FIGURE  12.   Plot  of  Terrain  Pressure  Versus  Forecast 
Minute  at  20N  HOE  (Experiment  IV) 
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FIGURE  13.   Plot  of  Terrain  Pressure  Versus  Forecast 
Minute  at  *»0N  HOE  (Experiment  IV) 
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FIGURE  Ik.      Plot  of  Terrain  Pressure  Versus  Forecast 
Minute  at  80M  HOE  (Experiment  IV) 
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FIGURE  15.   Plot  of  Terrain  Pressure  Versus  Forecast 
Minute  at  20N  110E  (Experiment  V) 
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FIGURE  16.   Plot  of  Terrain  Pressure  Versus  Forecast 
Minute  at  4 ON  HOE  (Experiment  V) 
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FIGURE    17 
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Plot  of  Terrain  Pressure  Versus  Forecast 
Minute  at  80N  HOE  (Experiment  V) 
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CHART  A.   Initial  Surface  Pressure  Analysis  From 
FNWC  Tape.  12  Z  10  May  73 
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CHART  B.   Surface  Pressure  Analysis  From  Model 
Interpolation.   12  Z  10  May  73 
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CHART  C.   12-Hour  Surface  Pressure  Forecast  Using 
a  =  0.0  and  Euler  Backward  Every  6  Hours 
VT  00  Z  11  May  73  (Experiment  I  &  II). 
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CHART  D.   12-Hour  Surface  Pressure  Forecast  Using 
a  =  0.3  and  Euler  Backward  Every  6  Hours 
VT  00  Z  11  May  73  (Experiment  I) 
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CHART  E.   12-Hour  Surface  Pressure  Forecast  Using 
a  =  0.4  and  Euler  Backward  Every  6  Hours 
VT  00  Z  11  May  73  (Experiment  I) 
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CHART  F.   12-Hour  Surface  Pressure  Forecast  Using 
a  =  0.5  and  Euler  Backward  Every  6  Hours 
VT  00  Z  11  May  73  (Experiment  I) 
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CHART  G.   12-Hour  Surface  Pressure  Forecast  Using 
a  =  0.0  and  Euler  Backward  Every  1  Hour. 
VT  00  Z  11  May  73  (Experiment  II) 
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CHART  H.   12-Hour  Surface  Pressure  Forecast  Using 

a  =  0.0  and  Euler  Backward  Every  2  Hours. 
VT  00  Z  11  May  73  (Experiment  II) 
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CHART  I.   12-Hour  Surface  Pressure  Forecast  Using 
a  =  0.0  and  Euler  Backward  Every  3  Hours 
VT  00  Z  11  May  73  (Experiment  II) 
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CHART  J.   12-Hour  Surface  Pressure  Forecast  Using 
Iteration  Initialization  (One  Hour) 
VT  00  Z  11  May  73  (Experiment  IV) 


60 


CHART  K.   12-Hour  Surface  Pressure  Forecast  Using 
Iteration  Initialization  (Three  Hours) 
VT  00  Z  11  May  73  (Experiment  IV) 
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CHART  L.   0.9  Sigma  Level  Initial  Guess  Stream  Function  Field 
(Experiment  V) 
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CHART  M.   925  MB  Level  Initial  Guess  Stream  Function  Field 
(Experiment  V) 
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CHART  N.   0.9  Sigma  Level  Final  Stream  Function  Field 
(Experiment  V) 


6^ 


CHART  0.   925  P^B  Level  Final  Stream  Function  Field 
(Experiment  V) 
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CHART  P.   0.9  Sigma  Level  Isotach  Field 
(Experiment  V) ' 
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CHART  Q.   925  MB  Level  Isotach  Field 
(Experiment  V) 
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CHART  R.   0.5  Sigma  Level  Initial  Guess  Stream 
Function  Field  (Experiment  V)  ■ 
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CHART  S.   500  MB  Level  Initial  Guess  Stream  Function  Field 
(Experiment  V)' [ 
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CHART  T.   0.5  Sigma  Level  Final  Stream  Function  Field 
(Experiment  V) 
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CHART  U.   500  MB  Level  Pinal  Stream  Function  Field 
(Experiment  V)!") 
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CHART  V.   0.5  Sigma  Level  Isotach  Field 
(Experiment  V) 
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CHART  W.   500  MB  Level  Isotach  Field  (Experiment  V) 
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CHART  X.   12-Hour  Surface  Pressure  Forecast  Without 

Terrain  Using  Sigma  Surface  Initialized  Winds 
(Experiment  V) 
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CHART  Y.   12-Hour  Surface  Pressure  Forecast  Without 

Terrain  Using  Pressure  Surface  Initialized  Winds 
(Experiment  V) 
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CHART  Z.   72-Hour  Surface  Pressure  Forecast  Without 

Terrain  Using  Sigma  Surface  Initialized  Winds 
(Experiment  V) 
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CHART  AA.   .9  Sigma  Level  Initial  Guess  Stream  Function 
Field  (Experiment  VI) 
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CHART  BB.   12-Hour  Surface  Pressure  Forecast  Without 

Terrain  Using  Sigma  Surface  Initialized  Winds 
VT  15  Z  10  Kay  73  (Experiment  VI) 
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CHART  CC.   72-Hour  Surface  Pressure  Forecast 

VT  12  Z  13  May  73  (Experiment  VIII) 
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CHART  DD.   FNWC  Objective  Surface  Pressure  Analysis 
12  Z  13  May  73  (Experiment  VIII) 
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CHART  EE.   84-Hour  Surface  Pressure  Forecast 

VT  00  Z  14  May  73  (Experiment  VIII) 
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CHART  FF.   FNWC  Objective  Surface  Pressure  Analysis 
00  Z  lH   May  73  (Experiment  VIII) 
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APPENDIX  A 

SIGMA  COORDINATE  -  LINEAR  BALANCE  EQUATION  IN 
FINITE  DIFFERENCE  FORM 


^  +  VjL^Vf  =  jrVJ,  +  Vir^Vj,  +  R   (TV2u  +  VT  .  V7t) 
(1)    (2)     (3)    (4)       (5)       (6) 


(1)       V2*  =  (*±J+1  -  2^  +  ty^J  -  Bj(*1J+1  -  ♦   x) 


(2)      —^-   -  ■#  (f-+1  -  f-)_l)(*iJ+l  "  ^ij-l5 


,2  .    lft 
s 


(3)      ffii.^[(#       2^  +^_x) 


(5)   -f  TV2,  =  -^  [(ulj+1  -  2w±J  +  ir^) 

"  Bj(7Tij+l  "  "ij-l5   +  AJ(7Tl+ij    "   27rij    +  TTi-ij)] 
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(6)         H  7T.V7r  .  ^R  [(T  .  T         )<-  -  ,         ) 


+  Aj(Ti+iJ   "  Ti-ij)(,,i+iJ   "  "i-lj'l 


A       =    ± 

1  2 

J      cos  e 


n      _  AG  tan   6 
Bj  2 


F.    =   2Q   sin   0 
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APPENDIX  B 
COMPUTER  PROGRAM 
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